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AXISYMMETRIC EXPANSION OF A PLASMA IN A MAGNETIC 


NOZZLE INCLUDING THERMAL CONDUCTION 
by Eddie L. Walker and George R. Seikel 
Lewis Research Center 

SUMMARY 

The axisymmetric expansion of a hot-electron, cold-ion plasma in the m^netic 
nozzle produced by a pair of Helmholtz coils is calculated with thermal conduction in- 
cluded. The plasma is assumed fully ionized (in the sense that neutrals have no signif- 
icant effect on the expansion) and the electron Hall parameter is assumed to have an ef- 
fective constant value. The analysis is based on the reduction of the partial differential 
kinetic equations to ordinary differential equations which are valid near the axis of sym- 
metry. The resulting nonlinear coupled equations were solved numerically by starting 
far downstream with the leading terms of an analytic asymptotic solution and working up- 
stream. The numerical results depend on two parameters: (1) the location of the on- 
axis sonic point, and (2) a dimensionless magnetic field /3, which involves the ratio of 
the maximum value of the applied magnetic field to the effective electron Hall parameter. 

The analysis predicts the spatial variation along the axis of symmetry of the ion 
flow velocity, electron temperature, plasma potential, and electron number density. 

Also presented is the dependence on ]3 of the spatial variation of the on-axis ion flow 
velocity and electron temperature. The analysis appears to be consistent with available 
experimental data. 

The final results give the on-axis ratio of the directed ion exhaust energy to the 
maximum (upstream) electron thermal energy as a function of the sonic point location 
and the dimensionless magnetic field The effect of including thermal conduction is 
striking; the energy ratio attains values more than an order of magnitude greater than 
that for an adiabatic expansion. 

The analysis is valid for a finite range of the parameter The largest value of 
/3 is limited by the stated assumption that the ion temperature is negligible, while the 
smallest value is limited by the implicit assumption of small electron cyclotron radii. 
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INTRODUCTION 

Those performing research on magnetoplasmadynamic (MPD) arc thrusters have 
for a number of years been investigating the steady-state acceleration of plasma in a 
magnetic nozzle. In an early paper (ref, 1) the following model of this process was de- 
duced for low-density dc discharge thrusters (i.e. , long ion mean free path); 

The electric power of the discharge is primarily added to the random 
energy of the plasma electrons. The energetic electrons ionize the 
propellant, and, as the plasma’s electron gas expands out of the de- 
vice, its random electron energy is converted to directed energy. 

Since there can be no divergence of current, however, the electrons 
drag the ions along. Physically this is accomplished by an axial 
electric field that the plasma itself builds. This field retards the 
electrons' expansion, accelerates the ions, and causes an azimuthal 
electron current to flow that provides the electromagnetic reaction 
force on the accelerator. The energy added to the ions is at the ex- 
pense of the random electron energy. The expansion process is con- 
trolled by the magnetic nozzle action of the spatially varying magnetic 
field. 

The existence of this plasma generated axial electric field was supported by the ex- 
periments of Domitz (ref. 2) and the research of Meyerand et al. (ref. 3). 

In a later paper (ref. 4) it was shown that on the axis of symmetry this axial elec- 
tric field was given by the gradient of the electron pressure divided by the electron 
charge density. Subsequent experiments by Bowditch (ref. 5) show that this relation for 
the axial electric field was applicable in the exhaust field of a fully ionized, low-power 
thruster. Hence, his results experimentally substantiated the proposed acceleration 
model. 

A detailed examination of the data of reference 5 shows, in addition, that the plasma 
expansion process is not adiabatic. The exhaust ion energy measured is much higher 
than one would obtain in an adiabatic expansion. And the electron temperature does not 
spatially decrease as rapidly as one would expect in an adiabatic expansion. 

For the electron temperature range of 5 to 10 electron volts measured by Bowditch, 
a plasma is a good thermal conductor. Thus, thermal conduction can be important. 

The object of this study is to analyze the expansion of a plasma in a magnetic nozzle 
including thermal conduction. Specifically, the expansion of a fully ionized plasma con- 
taining hot electrons and cold ions is examined in the magnetic nozzle produced by a 
Helmholtz set of coils, (The phrase ’’fully ionized” simply means here that any neu- 
trals which may be present play no appreciable role in the expansion.) The analysis is 
restricted to the vicinity of the axis of symmetry. The set of equations to be solved 
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includes the electron and ion continuity and momentum equations, Maxwell’s equations 
and an electron energy equation. As indicated in reference 6, the solutions can be used 
in an energy balance which leads to setting upper bounds on the potential performance of 
applicable MPD arc thrusters. Finally, the limits on applications of the analysis are 
discussed along with the results. 


THEORETICAL ANALYSIS 

In this section Maxwell's equations, the electron and ion continuity and momentum 
equations, and the electron energy equation are combined to obtain two coupled nonlinear 
ordinary differential equations for the electron temperature and ion flow velocity. The 
equations are valid for points near the axis of symmetry, with the effective electron 
Hall parameter and the magnetic field appearing as parameters. 


Assumptions 

The coordinate system for the plasma expansion is shown in figure 1. In order to 
reduce the usual kinetic equations to ordinary differential equations it is assumed that 
each variable of interest can be expanded in a power series in the radial coordinate r 
near the axis of symmetry. Then for instance, the electron temperature T^ is given by 

Te = Tg(r, z) = T^°^(z) + rT^^\z) + • (1 ) 

so that the electron temperature on the axis is simply 

Tg(r = 0,z) =T^°)(z) (2) 

By expanding each variable in this manner, substituting the series into the partial differ- 
ential kinetic equations, and comparing coefficients of like powers of r, ordinary differ- 
ential equations are obtained for the electron temperature and ion flow velocity on the 
axis. 

In deriving the final equations, the following additional assumptions are made (sym- 
bols are defined in appendix A and all equations and quantities are in S.I. units): 

(1) Steady state exists, 3/9t = 0. 

(2) Neutrals have no significant effect on the expansion (equivalent to assuming a 
fully ionized plasma). 


3 



(3) The species pressures are scalar; that is, viscous effects are negligible. 

(4) T. « Tg. 

(5) Quasi-charge neutrality exists on the axis, 

(6) There is no axial current density on the axis, = 0. 

(7) The on-axis induced magnetic field is negligible compared to the applied field. 
(Appendix B contains a discussion of the applied magnetic field.) 

(8) The electron inertia term in the electron momentum equation is neglected. 

(9) Inelastic collisions are negligible. 


Kinetic Equations 


The conduction current density is given by 

J=.q(N.u. -N^Ug) 

From the assumptions 



and 


it then follows that 


n(o) = n(0) 


u(0) = u(0) 

e,z i,z 


(3) 


(4) 


(5) 


( 6 ) 


Continuity equations. - The electron continuity equation is 


from which, 


u(») = 0 
e,r 


(8a) 


e,r 


_L ^ 

e 


e e,z 


(8b) 


where it has been assumed that 0. Completely analogous results are found for 

the ion flow velocity coefficients : 


« = 0 


(8c) 


u(l) . 

i,r 


i_± nP'uS") 

(0) dz L ^ 


2N>' 


(8d) 


Momentum equations. - The electron and ion momentum equations are, respectively. 


VP^H-qN^(EH.u^xB) = T^ 

(9) 

Nj^mj^u^ • Vuj^ - qNj^ + Uj^ x B ^ = - 1^ 

(10) 


where the ion scalar pressure and the gravitational acceleration have been neglected. 
Since Eg = 0 (see appendix C, eq. (CIS)), the 0- component of the sum of equations (9) 
and (10) is 


(“i ■ 0 - (-^z^r - "^r^z) = 0 

2 

Since is of order r and is of order r (see eqs. (CIO) and (C24), respec- 
tively), it follows that, to first order in r. 


(u.-Vu,)g=0 (12) 

2 

Now from equation (12) it is possible to show that u. ^ is of order r . First, using 
the assumption 9/90 =0 yields 
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( 13 ) 


(“i ‘ ’^i) e = “i r • +“i Z — ^ 

V 1 1/ ^ ^ 9r r ^ 3z 


9u, 


Then, the substitution of the r-expansions into equation (13) and the use of equation (12) 
give 


2„p),aP> = 0 

i,r 1 ,/. 


(14) 


where the fact that = 0 has been used. (The azimuthal component of velocity must 

vanish on the axis.) Substituting the result (8d) into equation (14) yields 


dlnuif) dl.rN(°)u(°)] 


dz 


dz 


(15) 


The solution of equation (15) is simply 


uS^KnPuS^ 


(16) 


where K is a constant. The ions are assumed to be monoenergetic, being created at 
some plane z = z* with no azimuthal velocity: 


Ui^g(z = z*)=0 


(17) 


It then follows from the expansion 


+r\(2)^(z) + 


that 


u[^0(z*) = 0 

and hence the constant K in equation (16) must vanish, assuming that 


Nfu(») 


z=z* 


* 0 
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Finally, 




( 18 ) 


The r-component of the Ion momentum equation (10) is 


N.m. u. 


1 i\ i,r 


3u. uf J 
— ii£ - + u, 

r 


9u, 


i,z gj y “i, 0®z " '^i, ~ ^e,r 


(19) 


which to first order in r becomes 




( 1 ) 

"i,r 


U: 


\u(0) ^ 

dz 


aN^O) M ^ j(l) 
qNi Ej, Ig^j, 


(20) 


where equations (8c), (18), and the results of appendix C, equations (C5) and (C20), have 
been used. Substituting the results (8d) and (C21) into equation (20) gives 


“i^i 


( 0 ) 


[n(0)] 


2 dz 


[n(»)u<0)]' 


U> 


(0) ^ j 1 d 

^ dz 1^(0) dz 


n(0)u<0) 

1 1, Z 


2 dz 


( 21 ) 


The z- component of the electron momentum equation (9) is 


dz 


+ qNe(Ez+Ue^rB0-Ue,0Bp=Ig^2 


( 22 ) 


which, to zero order in r, is 




dz 


(23) 
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The z-component of the sum of equations (9) and (10) is 


miNi 



3u. 


i,z 


9r 


+ a 


i,z 



9P^ 


( 24 ) 


which to zero order in r becomes 


du.(®^ 


dP 


,( 0 ) 


1 1 


dz 


dz 


= 0 


(25) 


Electron energy equation . - The final kinetic equation to be considered is the elec- 
tron energy equation 



VPe 


+ - ^ 
2 ® 


u. + V 




(26) 


where and 1^, are, respectively, the electron heat flow and the electron energy 
collision integral. To zero order in r, equation (26) becomes 


3 „(0) . 5 p(0) 

2 dz 2 ® 


L dz 




dQ 


( 0 ) 

e,z 


dz 


- I 


( 0 ) 

E 


(27) 


(As is shown in appendix D, the radial electron heat flow vanishes on the axis, so that 

V = 0.) Substituting equation (8b) into (27) gives 
e, r 


3 

2 


u 


( 0 ) 

e, z 


dP 


( 0 ) 


dz 


>(0)u(0) 


e,z 


d In N 
dz 


( 0 ) 


+ 2Q« * 


dQ 


( 0 ) 


e,z _ j(0) 

d~ E 


In order to close the set of equations (5), (6). (21), (23), (25), and (28), it remains 

to determine the collision integral coefficients and the electron heat 

/i\ e,r e,z e 

flow coefficients 

G f ^ G y r 


Transport Coefficients 

Momentum collision integral coefficients. - The diffusion Mach number, which is 
defined as 


8 



e = 


u. - u 
1 1 e' 


2k ( -1 + -£ 

V“i “ey 


nl/2 


is zero on the axis of symmetry by virtue of equations (6), (8a), (8c), and the fact that 
the azimuthal components of the electron and ion flow velocities must vanish at r = 0. 
From this fact it is reasonable to assume that £ is small in the vicinity of the axis of 
symmetry. Then, based on the well-known Grad-thirteen moment approximation for the 
electron and ion distribution functions, the electron momentum collision integral is, to 

first order in 

Qj/Pi|/2kTi/mi, 


€ and in the dimensionless heat flows /p.i/ik 


e' e 


and 


e e e ei 


(“1 - “ e ) 




(29) 


Here is the effective collision frequency for transfer of momentum between elec- 
trons and ions, and y is a dimensionless number of order unity involving the inter- 
particle force law (for the usual Coulomb force, y = 0.6). The derivation of equation 
(29) is a straightforward but tedious process; the details can be found in reference 7. 
The result depends in part on the assumption that « Qg/m^. 

Substituting the results (eqs. (5), (6), and (8a) to (8d)) into equation (29) yields 


,( 0 ) . 


e, z 


,( 0 ) 


(30a) 


^1° r = ° 
e,r 


(30b) 


j(l) ^”^e^ e ^ei ^e, r 

p(0) 


e, r 


(30c) 


The result (30b) depends on the fact that = 0, as shown in appendix D. 

6, r 

Energy collision integral coefficient . - From reference 7, the electron energy 
collision integral at r = 0 is 


I 



IH 



n<»)A 

e ei 



(31a) 


where note has been made of the fact that the diffusion Mach number e is zero on the 
axis of symmetry. Then, since expression (31a) becomes 




„(0)p(0) 

ei e 


(31b) 


Electron heat flow coefficients , - Based on the same assumptions made in the dis- 
cussion of the momentum collision integral coefficients, the z-component of the elec- 
tron heat flow to order zero in r is 


Q. 


( 0 ) _ 

e,.z 


.( 0 ) 


dkT 


( 0 ) 


e ei 


dz 


(32a) 


where is a dimensionless number involving the interparticle force law (for the 
Coulomb force, S 1. 86), The r- component coefficients of are 


Q! 


(0) _ 


= 0 


(32b) 


Q, 


( 1 ) _ 


[W'] 


e,r 


f ^ 


e ,.m(^) 


qB 


kT' 

(0) e 






»( 1 ) 


few 


/Jlr'T'(O) 

p(0) 

® dz 


(32c) 


The results (32a) to (32c) are derived in appendix D. Implicit in their derivation is the 
use of the "transport approximation" wherein the mean time between collisions for the 
electrons and the electron mean free path are assumed to be small compared to, respec- 
tively, the characteristic time and distance intervals for significant changes in the non- 
Maxwellian part of the electron distribution function. In expression (32c) the electron 
cyclotron frequency is given by 
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(33) 


Ct> 


_qB 


m 


e 


and 


^e=' 


ei 


(34) 


It is noted that equation (32c) contains two new unknowns, and Hence, in 


(2) ,„H ,,(1) 

e 

order to close the set of equations (5), (6), (21), (23), (25), and (28), expression (32c) 
will not be used. Instead, the following assumption is made: in any r-z plane the 
electron heat flow is parallel to the electron flow velocity 




or 


“e r 

Q = Q „ 

'^e,r ^e,z 

e,z 


(35) 


Then from equations (8a) and (35), it follows that 


Q? r = ® 
e, r 


(36a) 


Q 


(1) - *'e?r ^(0) 
e,r ( 0 ) e,z 
e,z 


(36b) 


The effective electron Hall parameter . - To complete the evaluation of the necessary 
collision integral and heat flow coefficients only requires the determination of the colli- 
sion frequency In preliminary calculations, it was found that use of the classical 

expression for u^) (see eq. (3. 94) of ref. 7) gave results which were in serious disa- 
greement with experimental data. As a result of this the classical expression for t'V' 
was replaced by an effective value. In particular, the following assumption was made: 
the electron Hall parameter can be approximated by an effective value which is constant, 
at least near the axis of symmetry, 
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(37) 


(jO 


( 0 )^( 0 ) ^ 

e e 


(t) 


( 0 ) 

e 

,( 0 )' 

ei 


constant 


where 


0 ,( 0 ) = <^(r = 0,^) , (31 

^ ”e “e 

This assumption is based on an analysis by Solbes (ref. 8) which predicts (using quasi- 
linear plane-wave analysis) that a saturation in the effective Hall parameter can result 
from electrothermal instabilities. In reference 8, Solbes also compares his analysis 
with MHD generator experimental data and obtains good agreement. The assumption 
(37) also permitted an analytic asymptotic solution (i.e. , large z); no such solution 
could be found when the classical expression for was used. 


Final Equations 

When expressions (32a) and (36b) are used, the momentum collision integral coeffi- 
cients (eqs. (30a) to (30c)) become, respectively. 


e,z 


= 

e,r 


SrN^O) 

(39a) 

2o^ ® dz 


r = 0 

e,r 

(39b) 

5ri.T(0)^^e %,r 

(39c) 

2 £ ^ dz (0) 


e,z 



Then, when equations (5), (6), (8b), (8d), (C21), (31b), (32a), (36b), (37), and (39a) to 
(39c) are used, equations (21), (23), (25), and (28) become, respectively. 


mJ^ — (Nu)] 

UN^Ldz J 



^^dE_ 5y^^e /dlnu ^ dlnN \ ( 4 ^) 
dz 2 j? dz \ dz dz / 
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J 


— (NkT ) + qNE =-lZ 
dz ® 2 ^ dz 


m.Nu ^ + A (NkTj = 0 
^ dz dz ® 


(41) 

(42) 


3 

d(kT ) 
u - 

MkT 

kTg d(kTg) 

d In u 

.A. 

/kT^ dkT A' 

2 

dz 

® dz 2 Jq 

B dz 

dz 

dz 

\B dz J 


where 


+ 3qB » 

m. (co^T ) 

1 e e' 


(43) 


N = = Np^ 

e 1 


„ , JO) . (0) 

e j z ^ 


T 

e e 


E = E<“> 


B = B(r = 0,z) 


c =S§ 
^ “e 


(44a) 

(44b) 

(44c) 

(44d) 

(44e) 

(44f) 


r, = T»)=J_ 

,.<0) 

ei 


(44g) 


Equations (40) to (43) constitute four independent equations in the four unknowns N, 
u, Tg, and E; the effective electron Hall parameter and the applied on-axis mag- 

netic field (which is approximately equal to B) are assumed to be specified. The elec- 
tric field E and its derivative dE/dz can be eliminated by combining equations (40) 
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and (41). The number density N, which occurs throughout the equations in the form 
d(ln N)/dz, can be eliminated by using equation (42). The final equations for the elec- 
tron temperature on the axis T^ and the ion flow velocity on the axis u are then 


m/ 5,,'^dVe)W^5^'^ j /dkT,^ 2 


kT 


2lj 


dz‘ 


2kTg 2 £/ kTg y dz 


, 5y/, d^u 


2£ 


kT, 


dz dz 


dz2 


+ m. 


3 ^ m.u 


kT, 


1 + 


mm 

2kT, 


^dz) 


(45) 


dln(kTg) a^i'^ du (^eV 

11 ■ " 4 ” — — — — 4 * — 


dz 


5 kT dz 
e 


dfqB 


d In u d In B ^^e\ 


dz 


dz 


dz 


dz 


d^kTg)' 

dz2 




qB 


= 0 


5 m.(a) T ) 
i' e e' 


(46) 


For both convenience and a reduction in the numerical calculations, equations (45) and 
(46) can be put into dimensionless form. It is convenient to nondimensionalize the ion 
velocity by the final exhaust velocity, the electron thermal energy by the final ion di- 
rected energy, and the distance by the radius a of the magnetic field coils (fig. 1). 
Thus letting 


-1 

(47a) 

U- - u 

(47b) 

u(z = OO) 



T = 



(47c) 
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c=-5z 

3 S 


(47d) 


^ 3t^ Q3, 

V 1 oo e e 


(47e) 


equations (45) and (46) take, respectively, the following dimensionless forms; 




- 2(U’)^ 


1+ ^ 
2 V T 


1 H-3U! 

^ 2 r - 


(48) 


2 

T” = — i + C.B 

5f ^ 


^uii + 6 u!u:\ +T* +^.iL 


T 5 T 


U B T 


(49) 


where all derivatives are with respect to the dimensionless distance x (e.g, , 

2 2 '07 

r” = (d r/dx )). Using assumption (7) and expressions (B2) and (B3) for the applied on- 
axis magnetic field, equation (49) becomes 




2J 


T 5 T 


U F T 


(50) 


where 


13 = 


_VT J 


max a 


8 (t^^T„) V m. u^ 

' e \ 1 « 


(51a) 


is a dimensionless magnetic field, and 


F = (l +x^) + [l + (1 +x)^] 


-3/2 


F» = - 3|x (l +x2)‘^^^ +(1 +x)[l + (1 +x)2_ 


■5/2 


(51b) 

(51c) 
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In expression (51a) Bjjjgx maximum value of the applied magnetic field, occur- 

ring on the axis of symmetry midway between the two field coils (fig. 1). 

Equations (48), (50), and (51) thus comprise the dimensionless set of nonlinear 
equations to be solved. 


NUMERICAL SOLUTIONS 

To the authors* knowledge equations (48), (50), and (51) do not admit anal 5 rtic solu- 
tions. However, it is possible to obtain asymptotic anal^ic solutions for large x by 
expanding t and U in power series in the variable x“ . These solutions represent 
relations between the starting values used for the numerical solutions. Then, instead 
of needing four initial values (i.e. , t, U, t', U'), it is only necessary to specify two 
initial values in order to begin the numerical solution at some point far downstream. 

In addition to deriving the asymptotic analytic solution in this section, the numeri- 
cal procedure is also discussed. Included in this discussion is the way in which experi- 
mental data were used to estimate representative values of the various parameters in 
the problem. This eliminated the potentially difficult problem of hunting numerically 
for values of these parameters that would lead to meaningful solutions. 


Asymptotic Analytic Solutions 


The analytic solutions for large x are based on the following series expansions: 


T 



X 



(52a) 


U = 1 




X 



(52b) 


(52c) 


There are no zero order terms in expansions (52a) and (52c) since T^ (and hence t) 
and N are assumed to vanish as x approaches infinity. The zero order term in (52b) 
is obvious since U(=u/u^) approaches unity as x approaches infinity. In general, the 
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Al_ 

term of a series such as equation (52a) could be assumed to be of the form 

-n-72 

Aj^x ; however, in appendix E it is shown that here Y 2 must be zero. 

Substituting expressions (52a) and (52b) into equation (48) and equating the sum of 
the lowest order terms in x” ^ to zero yield 

Aj +4BjAj + 3Bj = 0 (53) 

for which there are two solutions : 

Aj = -Bj (54a) 

Aj = -3Bj (54b) 


Only one of these solutions is admissible. To see this, reference is made to equation 
(42), which In dimensionless form is 

3NUU’ + Nt* + tN' = 0 (55) 


Then to lowest order in x”^, 

D^(3Bj +2Aj) = 0 (56) 

From the results (54a) and (54b), it follows that 

Dj = 0 (57) 

Equation (55) then yields to lowest order in x“^ 

3D2(B^+Aj)=0 (58) 

so that 

A^ = -Bj (59) 

provided that 02=5^ 0. (If D 2 = 0, both eqs. (54a) and (54b) are incompatible with eq. 
(55) and no solution exists.) 
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Substituting expansions (52a) and (52b) into equation (50) and equating the sum of the 
lowest order nonvanishing terms in x"^ to zero gives 


*.-T 


(60) 


Substituting equation (59) into (60) yields 




(61) 


Then, from equations (59) and (61) the asymptotic solutions (52a) and (52b) become, 
respectively. 


T = 





U = 1 



(62a) 


(62b) 


The higher order coefficients (Ag, Bg, etc.) could be determined in a similar manner in 
terms of and however, the algebra quickly becomes prohibitive. Moreover, it 
is anticipated that the truncated series (62a) and (62b) will be adequate for the purpose 
of starting the numerical solution. 

From equations (62a) and (62b), the first derivatives are, in truncated form. 


t' 



(62c) 


U’ = - 



(62d) 


These results assume, of course, that the asymptotic series (52a) and (52b) can be 
differentiated term by term. 
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Numerical Procedure 


Assuming that equations (62a) to (62d) are sufficiently accurate expressions for 
large x, the numerical solution can be started provided the coefficient is known 
O is treated as a parameter) . For a given ^ and an arbitrary value of Bj, equations 
(48) and (50) are solved numerically using equations (62a) to (62d) as starting values and 
working upstream (i.e., decreasing values of x). The quantities obtained are Xj^, the 
value of X for which the Mach number of the flow is unity, and 


1/t, 


max 



m^^u 



e,max 


(63) 


For each value of the parameter ^ two curves are obtained, one with Xj^ as a function 
of Bi and one with 1 /t „„ as a function of B,. From these sets of curves, two f am - 
ilies of curves are obtained showing (1) the dependence of for fixed 

values of / 3 , and (2) the dependence of on |3 for fixed values of Xj^. Finally, 

by repeating the latter set of calculations, but with the electron energy collision integral 
included, it is possible to determine an upper limit on ^ beyond which the analysis is 
no longer self-consistent. 

As indicated previously, the quantity (eq. (63)) is the primary goal of this 

analysis since it gives the ratio of the final (downstream) ion directed kinetic energy to 
the maximum (upstream) electron thermal energy. 

Starting calculations . - In order to begin the numerical solution some estimate of 
the magnitude of B^^ for a given |3 is needed. Such an estimate can be obtained from 


the measurements of reference 5. In that reference. 


a = 3.6x10”^ meter 

(64a) 

= 2.5x10"^ Tesla 

max 

(64b) 

m. = 6.68x10"^® kg 

(64c) 


The final ion flow velocity, u_^^, of reference 5 can be calculated from the z- component 
of the ion momentum equation which is simply equation (42) minus equation (41): 
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(65) 


Illll■llllllllllllllllllllllllllllll 


d(u^) ^ ^ dV ^ 5 y 

2 dz dz 2 dz 

where V is taken as the on-axis plasma potential of reference 5 (i. e. , E = E^®^ = 
- dV/dz). Integrating equation (65) gives 




kT 


( 66 ) 


At the calorimeter survey plane of reference 5, the following measurements were made: 


z = 0. 14 meter 

(67a) 

2 

m.u 

^ - 72 volts 

2q 

(67b) 

kT^ 

— - = 5.2 volts 

q 

(67c) 

V = 26 volts 

(67d) 


The operation of substituting equations (67b) to (67d) into equation (66) and taking 

y = 0. 6, £ = 1. 86 (these values are those for a fully ionized gas as given in appendix D) 

yields 


2 

m.u 

= 94 volts (68) 

2q 

where it has been assumed that = 0 inasmuch as the plasma potential of reference 5 
is measured relative to the surrounding tank walls of the experiment. If the ions were 
completely collision-free, the result (68) would be replaced by the simple sum of (67b) 
and (67d), or 98 volts; the last term on the right hand side of equation (66), which arises 
from the electron momentum collision integral, accounts for the difference. 

From the result (eq. (68)), 


u^ = 2. 12x10^ m/sec 


(69) 
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Then from equations (62b), (64a), (67a), (67b), and (68), a rough estimate of is 

-0.485 (70) 

It should be emphasized that (70) is only an estimate inasmuch as there is no assurance 
apriori that the calorimeter survey plane of reference 5 is sufficiently far downstream 
to justify the use of the truncated series (eq. (62b)). This estimate is adequate, how- 
ever, since all that is required here is some rough idea of where to begin the afore- 
mentioned numerical solution scheme, that is, what reasonable value of B^ is needed 
for a given ^ value. 

The appropriate value of |3 is found by substituting equations (64a) to (64c) and (69) 
into expression (51a), 


„ _ 5. 3xl0~^ 

Substituting equations (64a), (67a), (70), and (71) into equation (62a) gives 

T = 0.104 + 3.52x10 ^ 


(71) 


so that from equation (68), 


kT 
e 

q 


2 

3 



T = 6. 5 + 


0.22 


(cOeTg) 


(72) 


For > 1, the result (eq. (72)) gives the following range for the electron tempera- 

ture at the calorimeter survey plane: 

kT 

6. 5 volts ^ — ? < 6. 72 volts (73) 

q 

This range differs from the measured result (eq. (67c)) by at most 30 percent. Part of 
this difference is undoubtedly due to the fact that the calorimeter survey plane of refer- 
ence 5 is not sufficiently far downstream to justify the use of either of the truncated se- 
ries (eqs. (62a) and (62b)), Nevertheless, the result (eq. (73)) is sufficiently accurate 
to permit the use of equations (70) and (71) as a ’’starting” point for the numerical solu- 
tion scheme previously outlined. 
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As various values of and ^ are used some idea of the starting value of x is 
required. In order that the asymptotic starting solutions (eqs. (62a) to (62d)) be reason- 
ably accurate it is necessary that the starting value x = be large enough to make 
the lowest order terms dominant. To accomplish this the following constraints are im- 
posed on equations (62a) and (62b); 


-B. 

i<0.05 


X 


o 


that is, 


Xq>-20Bj=Wj (74a) 





that is, 


and 


that is. 



> - 20 A z. w 

0 — 


(74b) 


(74c) 


(Note that Bj must be negative.) The requirement that the magnitude of each higher 
order term in equations (62a) and (62b) be less than 5 percent of the corresponding 


22 



dominant term is of course arbitrary. However, such a requirement should certainly 
ensure that these asymptotic starting solutions are reasonably accurate. The starting 
value of X is then taken to be 


= maxCw^.Wg.Wg) 


(74d) 


The derivatives (eqs. (62c) and (62d)) are then automatically accurate expressions pro- 
vided that the asymptotic series can be differentiated term by term. 

Numerical integration . - Equations (48) and (50) were solved numerically on an IBM 
7094 electronic digital computer using a fourth order Runga-Kutta technique (see, e.g. , 
ref. 9). In order to enhance the accuracy of the results, a variable step size. Ax, was 

XI- 

used. For the computation the step size was 


(Ax) 


0.002(Ax)j^_j 

9 - 1 

m-Z m-1 


(75) 


with an initial step size of (Ax)q = -0.01. The constraint (eq. (75)), which tends to main- 
tain AU at approximately 0.2 percent of the final value of U, namely unity, is again 
arbitrary. It should, however, certainly keep the truncation error of the Runga-Kutta 
scheme at a reasonably small value while also reducing the computer run time. The 
constraint (eq. (75)) also allows a level of accuracy not possible with a fixed step size 
near the upstream singularity point (U— 0), where U begins to change very rapidly. 

As already mentioned, one of the quantities to be calculated is Xj^, the value of x 
where the Mach number of the flow is unity. For this purpose the Mach number must be 
computed at each step of the numerical solution. For simplicity, the Mach number of 
the flow is based on the adiabatic sound speed of the gas which is given by 



where P and p are, respectively, the scalar pressure and mass density (on the axis) 
of the gas. Taking into account the fact that m^^ » m^ and T^ » T. , expression 
(76) becomes 
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so that the Mach number of the flow on the axis of symmetry is 


M=^ = 



u 


which in terms of U and r is 


M = 


3U 

(5t)1/2 


(77) 


The approximate value of Xj^ is found by a linear interpolation in the computer program. 


RESULTS AND DISCUSSION 


Singularity at Maximum Electron Temperature 


The results of a typical computer run are shown in figures 2 and 3. Here /3 = 0. 1 
and the sonic point is located at Xj^ = 1.09 (the coordinate x = z/a is measured from 
the center plane of the downstream coil of fig. 1). At the point where the dimensionless 
electron temperature t reached a maximum value there was an overflow in the com- 
puter storage and the computation was automatically stopped. Such an overflow occurred 
on every run and indicates a singularity point in equations (48) and (50). It will be re- 
called that the electron heat flow on the axis is, from expressions (32a) and (32b), 


Q 


( 0 ) 

e,z 


5 

2 


,(0) 


e ei 


dkT^°^ 

e 

dz 


(78) 


Now the computer results show that the singularity point (z = z^) is the point at which T^ 
attains its maximum value; that is, t' goes through zero at x = x^^. Hence, 
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dz 


> 0 for z < z. 


dz 


= 0 for z = z. 


dz 


< 0 for z > Zi 


so that from expression (78) 


^e% ^ ® for z < 


= 0 for z = z, 


>0 for z > zj 


(79a) 

(79b) 

(79c) 


The results (eqs. (79a) to (79c)) imply that there is a point heat flow source located at 
z = Zj, the existence of which results in the singularity. (The assumption of zero on- 
axis current density rules out the possibility of any distributed heat flow source.) The 
accompanying computer overflow results from the fact that the singularity forces the 
term du/dz to become indefinitely large as z approaches from the downstream 
side. The details of the singularity are presented in appendix F, 

The fact that the numerical solution cannot be carried through the singularity point 
is not catastrophic since T^ occurs at the singularity. Hence, 

the quantity of interest, can indeed be calculated. 


Axial Variation of Plasma Quantities 

Figure 4 shows profiles of the normalized ion flow velocity U, the electron number 
density N/N^^ (calculated from eq. (55)), the applied magnetic field and 

the dimensionless electron temperature t and plasma potential AV calculated from 
the dimensionless form of equation (66) : 
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(80) 


AV s, 


V - v„ 



= 1+^T-U^ 

3J 


At the point where the electron temperature reaches a maximum (just upstream of the 
sonic point) the ion flow velocity is rapidly increasing, while the electron number 
density and plasma potential are falling off sharply. The figure also shows that there 
is still significant ion acceleration relatively far downstream (e.g. , at three coil radii 
downstream of the front coil, the ion flow velocity has attained less than 85 percent of 
its final value). 


Spatial Dependence of U and t on )S 

Figures 5 and 6 show the behavior of the spatial variation of U and t as the di- 
mensionless magnetic field /3 is varied over two orders of m^nitude, with the sonic 
point fixed at ~ 0. (This point corresponds to the center plane of the downstream 
coil of fig. 1). It can be seen from these figures that the on-axis ion flow velocity and 
electron temperature approach their final values over a shorter distance as the param- 
eter /3 is decreased. Figure 6 shows that t increases with /3. 


Ratio of Final Ion to Initial Electron Energy 

Figure 7 shows the ratio of the final ion to initial electron energy 1/t „ and the 

sonic point location plotted against the parameter for a given value of the non- 
dimensional magnetic field /3. By combining a set of such curves the final numerical 
results shown in figures 8 and 9(a) are obtained. Figure 8 shows the variation of the 
ratio of final ion to initial electron energy 1 /t „ as a function of sonic point location 

for various values of /3. Figure 9(a) shows this energy ratio as a function of j3 for 
various sonic point locations. 

Figure 8 illustrates that is a strong function of sonic point location if the 

sonic point lies between 0.75 and one coil radius downstream of the center plane of the 
downstream coil. Figures 8 and 9(a) both show that is a strong function of /3 

for j3 between 10” and 10” and sonic point locations less than one coil radius. The 
fact that is a decreasing function of ^ for a given sonic point location does 

not imply that 1 /t^„„ is a decreasing function of B (see eq. (51a)). This is be- 
cause the "constant” effective value of the electron Hall parameter may well 
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depend on (The parameter ( ot ^ is "constant” only in the sense that, for a 

in 0 0 

given system, d(cOgTg)/dx ~ 0 for all x.) 

The curves of figure 9(a) all approach a limiting value of 1 /t » 1.66 as ^ 

becomes large. This asymptote is simply the adiabatic solution of appendix G, that is, 
no electron heat flow and no electron energy loss term. The asymptote serves as a 
check on the accuracy of the numerical results; it also indicates the marked dependence 
of the results on the electron heat flow. 


Limitations of Analysis 

Figure 9(b) contains a set of curves similar to those shown in figure 9(a). However, 
in the calculation of the curves for figure 9(b) the electron energy collision integral 1^, 
has been included (this corresponds to including the first term on the right hand side of 
equation (50), which is of order jS^). A comparison of figures 9(a) and (b) shows that 
the two sets of results are in close agreement for small values of with a maximum 
difference of about 25 percent at |3 = 0.2. Beyond this point, however, the two sets of 
results begin to differ sharply; at ^ - 1 no solutions exist for Xjyj ~ 0.7 m the calcu- 
lations which include Ig, while solutions exist for the full range of Xj^ up to about 
^ = 10 in the calculations which omit 1^,. 

From these results it is apparent that the electron energy collision integral Ij, be- 
comes an important factor for /3 > 0.2. Now Ig represents the random kinetic energy 
lost by the electrons to the ions (note that 1^, is negative, eq. (31b)); this transfer of 
energy causes the ion temperature to rise. Then, since the energy loss term becomes 
important for /3 > 0.2, it follows that the ion temperature must also become appreciable 
at this point. Hov/ever, the ion temperature has been completely ignored in this analysis 
and in particular it has been omitted from the electron energy collision integral itself 
(see eqs, (31a) and (31b)). Clearly then, the analysis is inconsistent for jS > 0.2. (The 
inclusion of I„ along with the simultaneous omission of T. results in a violation of the 
conservation of energy.) This fact, along with rough estimates of typical values of j3 
for low and high- power magnetoplasm adynamic arcs (e.g. , eq. (71)), suggests that the 
analysis herein is applicable to the former but not for the latter case. The ion temper- 
ature is expected to be important for the high- power arc. 

There is also a limit on the analysis for very small values of |8 which results from 
the implicit assumption that electron cyclotron radii are small compared to a character- 
istic length of the system (e.g. , the magnetic field coil radius a). This requirement is 
necessary in order for the electrons to cyclotron within the device. The limit can be 
illustrated by noting that from equations (51a) and (63) 
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1 


a 


(81) 


/3 



where a characteristic electron cyclotron radius has been defined: 


r 


ce 


8 ^^e, max 

177 m_ 


, 1/2 


^^max \ 
™e / 


(82) 


The numerator of expression (82) is simply the average speed for a Maxwellian velocity 
distribution with temperature T (ref. 10), while the denominator is the electron 

(3 ^ XXx 

cyclotron frequency at (Although usually does not occur at B eq. 

(82) nevertheless yields a typical value of the electron cyclotron radius.) If the ratio of 
the characteristic electron cyclotron radius to the coil radius is to be small, r^^/a « 1, 
then from equation (81), 


/3» 


/A£ 

V Ott 2 ym|^ 


From figure 8 the largest possible value of 
(83) becomes 



1 

y. 

max / 

is less than 32, 


(83) 


SO that expression 


|3 » 8. 5x10"^ 



V i / 


(84) 


The inequality (84) gives an indication of the lower bound on the dimensionless magnetic 
field /3. For the case in reference 5 (argon propellant), the requirement (84) becomes 




3xl0~^ 


( 85 ) 
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which is certainly satisfied by the experiment (see eq. (71)). 


Significance of Results 

The inclusion of thermal conduction in the expansion of a plasma in a magnetic 
nozzle yields a dramatic effect. As shown in figures 8 and 9(a), the ratio of final ion to 
initial electron energy can be increased more than an order of magnitude over its adia- 
batic expansion value. 

The results of this analysis should permit one to predict the final exhaust velocity 
and spatial variations of the plasma quantities of a low power MPD arc thruster as a 
function of upstream electron temperature. However, to do so requires that one know 
both the location of the sonic point in the flow Xj^ and the value of the dimensionless 
magnetic field j3. Although reasonable estimates of Xjyj and ^ can be made for a given 
device, there is at present no precise method of calculating these quantities. Thus, to 
compare the analysis herein with experiments requires very detailed diagnostics of the 
experimental plasma exhaust. The only experiment of this type is that by Bowditch (ref. 
5). Unfortunately, Bowditch was unable to obtain data in the neighborhood of the sonic 
point because the probes were destroyed by the large heat transfer in that region. Hence, 
presently, no quantitative comparison can be made between the analysis herein and ex- 
periment, There does appear, however, to be qualitative agreement. Using Bowditch's 
data to begin the numerical solution resulted in rapid convergence of the calculations; 
furthermore, the asymptotic analytic solution for electron temperature (eq. (73)) agreed 
relatively well with Bowditch's measured value. Thus, the analysis at least appears to 
be consistent with available data. 

One important application of the results of this analysis has been indicated. The 
analysis leads to specifying the maximum range of the ratio of final ion to initial electron 
energy, Seikel et al. (ref. 6), indicate how this information can be used in a power bal- 
ance or a low power dc MPD arc to set an upper bound on both the possible efficiency 
as a function of specific impulse and the maximum specific impulse attainable. 


CONCLUDING REMARKS 

The axisymmetric expansion of a thermally conductive plasma in a magnetic nozzle 
has been analyzed. As an important step in the determination of an upper bound on the 
potential performance of applicable MPD arc thrusters, the on-axis ratio of the directed 
exhaust energy to the maximum (upstream) thermal energy was calculated as a function 
of two parameters: (1) the location of the on-axis sonic point, and (2) a dimensionless 
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magnetic field /3, which involves the ratio of the maximum value of the applied magnetic 
field to an effective electron Hall parameter. The results indicate that thermal conduc- 
tion can be the dominant process in the expansion inasmuch as the energy ratio can be 
significantly higher than that of an adiabatic expansion. The analysis appears to be valid 
for the low power MPD arc; the point where it becomes inconsistent indicates that the 
analysis must be extended by including the ion temperature to cover the case of the high 
power MPD arc (large 0 ). The implicit assumption of small electron cyclotron radii 
sets a lower bound on j3. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 16, 1970, 

120-26. 
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APPENDIX A 


SYMBOLS 


m 

th, 

coefficient in series for 
dimensionless electron 

k 

Boltzmann constant, 
1. 3805x10"^^ J/K 


temperature 

a 

length defined by fig. 10, m 

a 

effective radius of magnetic 
field coil, m 

M 

Mach number of flow 
electron rest mass. 


adiabatic sound speed of gas, 
m/sec 


9. 109X10"^^ kg 


™i 

ion rest mass, kg 

B 

magnetic flux density, T 
m"^ coefficient in series for 

N 

“3 

number density, m 

m 

normalized ion flow vel- 

0 

”of the order of*’ 


ocity 

P 

p 

scalar pressure, N/m 

C 

dimensionless force law par- 


2 

heat flow, J/m -sec 


ameter defined by eq. 
(47d) 

q 

electron charge, 
1. 602X10"^® C 

Cl 

parameter defined by eq. 
(47e), T"^ 

r 

radial coordinate, m 


random electron velocity, 

r 

unit vector in r-direction 

t; 

m/sec 

T 

temperature, K 

m 

m^^ coefficient in series for 

t 

time, sec 

- 3 

number density, m 

U 

normalized ion flow velocity 

E 

electric field intensity, V/m 

u 

flow velocity, m/sec 

F 

dimensionless function de- 


final ion flow velocity, m/sec 


fined by eq. (51b) 

V 

plasma potential, V 

% 

electron energy collision in- 

V 

total particle velocity, m/sec 

tegral, N/in^-sec 

Wi,W2,W3 

defined by eq. (74) 

h 

electron momentum collision 

X 

dimensionless axial distance 

T 

integral, N/m^ 
conduction current density, 

z 

axial coordinate, m 

K 

A/m^ 

constant defined by eq. (16) 

z’ 

dimensionless force law 
parameter 
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dimensionless magnetic field 
defined by eq. (51a) 

y 

dimensionless force law 
parameter 

Ax 

step size for numerical anal- 
ysis 

e 

diffusion Mach number 

^o 

permittivity of free space, 
F/m 

e 

azimuthal coordinate, radians 

9 

unit vector in 0- direction 


permeability of free space, 
H/m 

^ei 

effective electron- ion momen- 
tum transfer collison fre- 
quency, sec”^ 

£ 

dimensionless force law 
parameter 

P 

-3 

mass density, kg/m 

T 

dimensionless electron tem- 
perature 

^e 

V-ei 

OJ 

e 

Subscripts : 

electron cyclotron frequency, 
rad/sec 

e 

electron 

ext 

applied field 

i 

ion 

ind 

induced field 

M 

Xjyj is value of x where Mach 
number is unity 


m step number in numerical 

analysis 

max maximum value 

r radial component 

z axial component 

o starting value for numerical 

solution 

6 azimuthal component 

Superscripts: 

(m) m*' -order term with respect 

to r 

* z* is plane where ions are 

created 
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APPENDIX B 

THE APPLIED MAGNETIC FIELD 

A typical magnetic field coil configuration of MPD arcs is shown in figure 10. It is 
assumed here that the two field coils are closely wrapped and sufficiently narrow so that 
the magnetic flux density near the axis can be accurately determined on the basis of two 
equivalent single turn loops. Each such loop carries a current equal to the current in 
the actual coil multiplied by the number of turns of the actual coil, and has a radius (the 
dimension a of fig. 10) such that the magnetic flux density of the two loops is approxi- 
mately equal to that of the actual coils, at least near the centerline. 

For simplicity, it is assumed that a = f in figure 10, as is the actual case in ref- 
erence 5. The magnetic flux density on the centerline due to the external coils, 
is then (see ref. 11) 



where B^„^ is the maximum value of the applied magnetic field on the axis, occurrir^ 

XIIaX 

midway between the coils, at z = - (f/2) = - (a/2). Substituting x = z/a into express- 
ion (Bl) yields 

Bex.« = ^ B^ax ^ '"i <B2) 

and then 

^^extff ^ 5^ 3 |_ 3 ^ (j, ^^ 2 ) ^ ^x) [1 + (1 +xf] (B3) 

dx 16 L J 

The anal 3 dic asymptotic solution of equation (50) requires the evaluation of the ex- 
pressions in braces in equations (B2) and (B3) for large x. The binomial expansions of 

F(x) = \l+x^) + [l+(l+x)'^J (B4) 
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and 


= - 3x(l +x^) - 3(1 +x) [l +(1 +x)2] 

dx 


-5/2 


are, respectively, 


F(x) = 2x“^ - 3x‘^ +3x"^ - ix"® +0 (x"'^) 

2 


and 


= - 6x"^ + 12x"® - 15x"® +0 (x"'^) 
dx 


Finally, from the symmetry of the configuration shown in figure 10 it 


that 


®0,ext ■ ° 


(B5) 


(B6) 


(B7) 

is obvious 

(B8) 
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APPENDIX C 


MAGNITUDE OF FIELD QUANTITIES 


Magnetic Field 


The induced magnetic field is governed by the equation 


the z- component of which 




r 0r 


(Cl) 


(C2) 


"'since 3/ 36 = 0, Integrating equation (C2) with respect to r yields 


®0,ind 


(?)/ 


rJ^(z,r)dr + 


h(z) 


(C3) 


The arbitrary function h(z) in expression (C3) must vanish identically in order that 
ind finite at r = 0. Hence, since = 0(r), assumption (6), it follows that 

B,_j„^=0(r2) (C4) 

The 0-component of the total magnetic field is then 


Be “Bj, ) 

since = 0 by symmetry (see appendix B). 

The divergence of the total magnetic field is 


(C5) 


V . B = 0 


(C6) 


from which 


1 a 

lJ.(rBj+. ^ 


r 0r 


dz 


(C7) 
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since d/dd=Oc Integrating equation (C7) with respect to r yields 


B 


r 



dB^{r,z) 

dz 


dr 

r 


(C8) 


where the arbitrary function g(z) must be identically zero in order that be finite on 
the axis. Hence, to first order in r, 


B 


r 


/r\ ^Bz(0,z) 

[ 2 / dz 


(C9) 


so that 



(CIO) 


and 


r( 1) ^ 1 dB(0,z) 

^ 2 dz 


(Cll) 


Electric Field 

Because of the assumption of steady-state conditions the curl of the total electric 
field is given by 


V X E = 0 (C12) 

Hence, E can be expressed in terms of the gradient of some scalar function 


E = -V$ 


the 0- component of which 


E 


8 ~ ~ 


r do 


(C13) 


(C14) 


36 



Then since d/dd = 0, 




= 0 


The divergence of the total electric field is 


V 


E = . N^) 


(C15) 


(C16) 


from which 


1 1 (rE J + !^ = (N. - N ) (C17) 

rar r 3z \e^j 

since 3/00 = 0. Integrating equation (C17) with respect to r, and recognizing that 
must be finite at r = 0, yields 

= - (;) /' " (v) 

Hence, to first order in r, 


Er = - 


dEjO,z) 

dz 


since assumption (5)o Then 


e(0) = 0 


and 


^ 2 dz 


(C19) 


(C20) 


(C21) 
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Conduction Current Density 


With the assumption of steady- state conditions the conservation of charge equation 
becomes 


V . J = 0 


(C22) 


from which 


l-i(rJ)+— ^ = 0 (C23) 

r dr dz 

since 9/96^0. Substituting the r- expansions for J and J into equation (C2 3), and 
using the assumption J; ^ = 0, yields 

= 4^^ = 0 (C24) 
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APPENDIX D 


THE ELECTRON HEAT FLOW 


The material in this appendix is based on reference 7. Equation (5.4) of that refer- 
ence is the equation for the species’ heat flow in a general gas system. For the case of 
electron heat flow in a fully ionized gas, the equation becomes 


(b XQ ) VT = 6(i m c^c \ - 1^6 (m e ) 

B ^ 2 m^ ® \2 ® ® ®/ 2 m„ V e e/ 


(Dl) 


where ^ = v - TT^, with v being the total particle velocity, and where 6 [(l/2)mgCgCg 
and 5(mgCg) are, respectively, the total electron heat flow collision integral and the 
total electron momentum- collision integral for elastic collisions (the 6-notation is that 
used in ref. 7). In equation (Dl) the traceless pressure term has been ignored, in accord 
ance with assumption (3). Equation (Dl) is based on the ’’transport approximation” to 
the Boltzmann equation. This approximation is applicable to plasmas which satisfy two 
conditions: (1) the species’ distribution functions are close to their Maxwellian forms, 
and (2) the small non-Maxwellian parts are slowly varying functions of the time and spa- 
tial coordinates in the sense that their characteristic times and lengths are much larger 
than the corresponding collision times or mean free path lengths. Condition (1) has al- 
ready been assumed in the calculation of the collision integral coefficients in the text. 
Insofar as equation (Dl) is concerned, and in view of assumptions (1) and (3), condition 
(2) can be replaced by the requirement that the characteristic length for the variation of 
the electron heat flow, be large compared to the electron mean free path length. 

The derivation of the collision integrals in equation (Dl) is extremely tedious and 
only the results will be given here. To first order in the diffusion Mach number, e, and in 

Q Q; 

the dimensionless heat flows and — — , the electron mom- 

Pg(2kTg/m^)^/2 P.(2kT/m.)^/2 

entum collision integral is (eq. (3.42) of ref. 7), 


6(mgCg) = m^Ngi^g. 


(Ui - Ug) + 


yQ^ 


(D2) 


where the assumption « Q^/m^ has been used. To the same level of accuracy 

the electron- electron heat flow collision integral is (eq. (3. 50) of ref. 7) 
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(D3) 


5' j”e‘=IS)ee 



where a term involving the "Coulomb logarithm" has been omitted, with the condition 
that In A_. » 1/2. Typical values of In A which depends on T and N , are 

CX C li C tZ/ 

about 10 (see ref. 12). The electron- ion heat flow collision integral is (eq. (3.45) of 
ref. 7) 


.2- 



5(1 - y)Pg (u. - Ug) + (7z* - 5y - 2)Q^ 


(D4) 


where z' is-a dimensionless number which is dependent on the interparticle force law 
(the y used here is the z of ref. 7) . The total electron heat flow collision integral is 
given by the sum of equations (D3) and (D4). Then from equations (D2) to (D4) the right- 
hand side of equation (Dl) is 


I = 


V . 

ei 


5yPe(u. -u^) +2£Q^ 


(D5) 


where 


£ ^5y + 1+^}^ - I z» 


(D6) 


The z-component of equation (Dl) is, to zero order in r. 


5 e ^ _ j(0) 


2 m_ dz 


- 

ei ^e, z 


(D7) 


where equations (C5), (CIO), and (6) have been used. Hence, from equation (D7), 


Q! 


(0) _ 


5 Pf 


(D8) 


The 6- component of equation (Dl) yields 
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B 


(D9) 


B ®r 

% r = --^h % z 

ct) B ^ B 
e z z 


since 3/90 = 0, while the r- component yields 


Q 


e,6 


B 


9kT^\ B. 

I . 5_£_£ 


w B, \ " 2 m_ 3r / B„ 

© z 


Substituting equation (DIO) into (D5) gives 


l0 = - 


fu . 

ei 


5yp 4 -=!£ 


P„ 3kT \ 2£Ba 


Substituting from equation (D5) into (Dll) gives 




£u . 


ei) 


\ 2 B 


*^ei 


\2 


Q. 


e B, 


5rPg (u. - U J 


e,z 


- afy^ (u. - U y 


„ P^ 3kT^ Bn 

-5J ^ ® i + 2J — Q„ , 

B, m^cj 9r 
zee z 


Finally, substituting equation (D12) into (D9) yields 


Q, 


e,r 


‘eM/B 


2 /\B 


- 1-1 




■ 5£ 


Pe B 




+ 5yP^ 




CO^T B, 

e e z 




Now, u^^^n = 0 in order for the limit lim u„ n 0 to exist; similarly, T^‘ 
e, o j._^Q e, (7 e 


( 1 ) _ 


(DIO) 


(Dll) 


(D12) 



0 in order 
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no loss of generality is incurred by using the asymptotic expansion (eq. (52b)). Simi- 
larly, let 

OO 

T=x (E7) 

j=0 


N=x^^^D.x"j (E8) 

j=0 


where 




>'3 > 0, ^ 0 (ElO) 


Note from expression (ElO) that is taken to be positive since N— 0 as x— «>, which 
is obvious from physical considerations. The result t— 0 as x— which is not as 
obvious, can be shown by reference to equation (55): 


3NUU’ + Nr’ + tN’ = 0 (Ell) 

If T does contain a zero order term, then clearly ^2 = 0* Then substituting equations 
(E6) to (E8) into (Ell) yields, to lowest order in x“^. 


(terms of order 




0 


(E12) 


where 


m > ^3 +2 (E13) 

Clearly, since the last term in equation (E12) is the dominant one, at least one of the 
factors y^, A^, must vanish, but this is not possible in view of equations (E9) and 
(ElO). Hence, a contradiction has been reached and it follows that t cannot contain a 
zero order term. The expansion (E7) then becomes 
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j=0 

where now 


72 > 0, Aq ^ 0 


(E14) 


(El 5) 


To determine r2> the expansions (E6) and (E14) are substituted into equation (48) 
which yields, to lowest order in x"^, 


72(72 + 2)A^x 


2 - 272-2 


+ 6n(n + l)A^Bj^x 


-72-n-2 


+ 9n^B^x"^”"2 = 0 
n 


n>l (E16) 

where it has been assumed that B is the first nonzero coefficient beyond B in the 
expansion for U. If 72 < n, then clearly the first term in equation (El6) is the lowest 
order term, which requires 


72 (r2 + 2)A^ = 0 

which is not possible in view of the conditions in equation (El 5). Next, if yg ^ then 
the last term in equation (E16) is the lowest order term, which requires 

B = 0 
n 

but this case must also be ruled out since was stipulated to be nonzero. The re- 
maining case is Y 2 = n which makes all the terms in equation (E16) of the same order, 
and yields 


(n + 2 )Aq + 6(n + l)A^^Bjj + 9nB^ = 0 


(E17) 


The two solutions of equation (E17) are 



B„, - 3B. 


n 


(E18) 
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It is not possible to explicitly determine the value of n inasmuch as equations (50) and 
(55) provide no information in this regard. However, by comparing the calculated val- 
ues of kTg/q with the measured values of reference 5, the '’choice” Bj # 0 can be 
shown to be the best one. Assuming the first solution in equation (E18) and considering 
only the lowest order x-dependent terms, the expansions (E6) and (E14) become, re- 
spectively. 


U « 1 + Bj^x"“ 


T 


« A^x 
o 


-n 





-n 


so that 


T ^ 



(1 - u) 


and 


kT 
e 

q 


2 

3 



(E19) 


Substituting equations (67b) and (68) into (E19) gives 


_e «23.4/-^^— \ volts (E20) 

q \n + 2/ 


at the survey plane of reference 5. Then, for instance. 


kT 

— ^ ^ lo8 volts for n = 1 

q 


kT^ 

— - ~ 11. 7 volts for n = 2 

q 

Since n/(n +2) is a strictly increasing function of n, it is clear that the calculated 
values of kT /q increase with increasing values of n. The measured value of kT /q 

G G 

is 
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— - = 5. 2 volts 

q 

so that clearly the choice n = 1 is the most suitable one. Note that the second solution 
in equation (E18) is not allowable since equation (Ell) is, to lowest order in x“^ using 
equations (E6), (E8), and (E14), 

®o [^"®n + + ^3)^0] "" ° 

which becomes upon substitution of the second solution in equation (E18), 

® 

A contradiction has thus been reached since yg, and are all nonzero. Hence, 
the expansion (E14) becomes, with yg = n = 1, 

T A^x" A^ = - 0 (E22) 

j=0 

which has the same form as the series assumed in equation (52a). 

Finally, substituting the first solution in equation (E18) into (E21) gives 

yg = 2 (E23) 

Note that the result (eq. (E23)) is independent of both n and B^. Hence, the expansion 
(E8) becomes 

00 . 

j=0 

which is identical in form to the series (eq. (52c)) inasmuch as the first coefficient in 
that series is shown to be zero (eq. (57)). 
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APPENDIX F 


THE SINGULARITY POINT 

As mentioned in the text, the singularity point can be traced to the fact that the term 
du/dz becomes indefinitely large as z approaches z^ from the downstream side (i.e. , 
as • z — Zj), where Zj is the coordinate where the electron temperature attains its 
maximum value. To see this, reference is made to equations (50) and (55). Noting that 

lim T* = 0 (through negative values) (Fl) 

X— Xj 


and 


lim T*' < 0 

X— Xj 


(F2) 


equation (50) becomes, in the limit as x— Xj, 


r“(xi) |3V(xi) 
^ 2 £ ^ 


9^F(x,) 


max X— X 


lim (U U’) + lim 
+ 


lim (i:^\ 


(F3) 


where it has been assumed that U remains finite as x 
then from expressions (F2) and (F3) it follows that 


lim 

X— Xj 



<0 


Xj. Now, if lim U' >0, 

X— x^ 


(F4) 


Then, from expressions (Fl) and (F4) it follows that 


lim 

X— Xj 



= +CO 


(F5) 
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Hence, either 


lim U’ = +» (F6) 

X— X j 


or 


lim U = 0 (through positive values) (F7) 


or both (F6) and (F7) occur. However, it can be shown from equation (55) that (F7) can 
only occur if (F6) also occurs. Equation (55) is, in the limit as x xj", 


Slim (OTI')+V 3 ^ 

X— Xj 


lim 

X— Xj 



0 


(F8) 


Hence, if lim U = 0, it follows that 

X— Xj 


lim U' = +00 

X— Xj 


since lim N' < 0. Hence, in any event, 

X— Xj 


lim U' = +« (F9) 

X— Xj 

The fact that U' increases indefinitely as x — Xj results in an overflow in the 
computer storage. Several of the computed results show that U actually goes through 
zero and U' becomes prohibitively large as t' approaches zero through negative val- 
ues. The fact that U becomes slightly negative, in contradiction to (F7), is simply due 
to the inaccuracy of the Runga-Kutta scheme at the singularity point with respect to the 
rapidly varying function U. 
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APPENDIX G 


THE ADIABATIC SOLUTION 


In the adiabatic expansion there is no heat flow and no energy loss term, 
these conditions equation (28) becomes 


3 

2 



2 ® dz 


= 0 


where 



Using the relation P^ = NkT^, equation (Gl) yields 


3 d In kT^ ^ d In N ^ q 
2 dz dz 

Equation (25), which remains unchanged, 

m.Nu^+^ = 0 
^ dz dz 

can be written as 

d In N du 

dz dz kT dz 

e 


Combining equations (G2) and (G4) yields 


5 

2 


dkT^ 
e 

dz 


+ m^u 


du 

dz 


0 


which in nondimensional form is 


T’ + - UU’ = 0 
5 


Under 


(Gl) 


(G2) 


(G3) 


(G4) 


(G5) 


(G6) 
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The solution of equation (G6) is simply 


T +— = constant = t „ +— U? 

5 5 ^ 


(G7) 


where Uj is the value of U at the point where t reaches its maximum value. 
Now, as X— «», T-^ and U— 1, so that from equation (G7) , 


T 3 

max 


(l-U?) 


(G8) 


From equation (G6) it follows that 


Uj = U(xj) = 0 


(G9) 


or 


U’(Xi) = 0 


(GIO) 


the latter implying U is a minimum at x., or both (G9) and (GIO) occur. In any event, 

O 


it is permissible to neglect Uj in equation (G8), so that, finally, 

1 


/ adiabatic 


5 

3 


(Gll) 


The results shown in figure 9(a) were obtained using equation (50) with the first term 
on the right side omitted, that is, no energy loss term. 


r‘» = /3F — (t» +- UU’) + T’ fe' +^- — \ 

2 t\5/ \UFt/ 


(G12) 


As /3 — oOj it follows from equation (G12) that 


T* +- UU’ - 0 (for all x) 
5 


(G13) 
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in order for the term containing j3 to remain finite. The consequence (eq. (G13)) is 
simply the adiabatic equation (G6), of course, and it follows that the asymptote of figure 
9(a) is simply the adiabatic solution given by equation (Gil). 


52 



REFERENCES 


1. Seikel, G. R. : Generation of Thrust - Electromagnetic Thrustors. In Electric 

Propulsion for Spacecraft, NASA SP- 22, 1962, pp. 19-24. 

2. Domitz, S.: Experimental Evaluation of a Direct-Current Low-Pressure Plasma 

Source. NASA TN D-1659, 1963. 

3. Meyerand, R. G. , Jr.; Salz, F. ; Lary, E. C.; andWalch, A. P. : Electrostatic 

Potential Gradients in a Nonthermal Plasma. In Proceedings of the Fifth Int. Conf . 
Ionization Phenomena in Gases, vol. 1, 1961, H. Maecker, ed. , North- Holland 
Publ. Co., 1962, pp. 333-342. 

4. Seikel, G. R. ; Bowditch, D. N. ; and Domitz, S. : Application of Magnet ic- 

Expansion Plasma Thrustors to Satellite Station Keeping and Attitude Control 
Missions. Paper no. 64-677, AIAA, Aug. 1964. 

5. Bowditch, D. N. : Investigation of the Discharge and Exhaust Beam of a Small Arc 

Plasma Thrustor. Paper .no. 66-195, AIAA, Mar. 1966. 

6. Seikel, G. R. ; Connolly, D. J. ; Michels, C. J. ; Richley, E. A.; Smith, J. M. ; 

and Sovie, R. J. : Plasma Physics of Electric Rockets. Proceedings of the Con- 
ference on Plasmas and Magnetic Fields in Propulsion and Power Research, NASA 
SP-226, 1969, pp. 1-64. 

7. Walker, E. L. : Transport Phenomena of General Non- Equilibrium Gas Systems. 

Ph.D. Thesis, Case Institute of Technology, 1967. 

8. Solbes, A. : Quasi-Linear Plane Wave Study of Electrothermal Instabilities. In 

Electricity from MHD, vol. I, Proceedings of a Symposium on Magnetohydrody- 
namic Electrical Power Generation, International Atomic Energy Agency, 1968, 
pp. 499-518. 

9. McCracken, D. D. ; and Dorn, W. S. : Numerical Methods and Fortran Program- 

ing. John Wiley & Sons , Inc., 1968, p. 325. 

10. Delcroix, J. L. (Melville Clark, Jr.; David J. BenDaniel; and Judith M. 

BenDaniel, trans.): Introduction to the Theory of Ionized Gases. Interscience 
Publishers, Inc., 1964, p. 30. 

11. Plonsey, R. ; and Collin, R. E. : Principles and Applications of Electromagnetic 

Fields. McGraw-Hill Book Company, Inc. , 1961, p. 208. 

12. Spitzer, L. ; Physics of Fully Ionized Gases. Seconded., Interscience Publishers, 

1967, p. 128. 


53 



Dimensionless electron temperature, 


® Magnetic field coMs 0 


Figure 1. - Plasma expansion in a magnetic nozzle. 



Figure 2. - Dimensionless on axis electron temperature t as function of dimensionless 
distance x for aimensionless magnetic field p= 10 ’^ and sonic point location 
= 1.09, 



Figure 3. - Normalized on axis ion flow velocity U as function of dimensionless distance 
x for dimensionless magnetic field p = 10'^ and sonic point location X|y^ = 1.09. 


Normalized Ion flow velocity, U = u/u, 





1 2 3 4 5 6 7 8 9 10 


Dimensionless distance, x = z/a 

Figure 4. - Axial profiles of plasma properties for dimensionless magnetic field 
3=^ 10'^ and sonic point location x^ = 1.09. 



0 123456789 10 

Dimensionless distance, x =* z/a 


Rgure 5. - Normalized on axis ion flow velocity U as function of dimensionless distance 
x with dimensionless magnetic field p as parameter. Sonic point location -0, 
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Dimensionless electron temperature, 



Figure 6. - Dimensionless on axis electron temperature t as function of dimension- 
less distance X with dimensionless magnetic field p as parameter. Sonic point 
location, x^y| = 0. 



Figure 7. - Ratio of final ion to initial electron 
energy Ixation x^ 

as function of parameter Bj for dimension- 
less magnetic field p = 5xl0"2. 


Asymptote 



.xatio 01 final ion to initial electron energy function of sonic point location with dimensionless magnetic field [J as parameie- 






z = 0 z 

Figure 10. - Typical configuration of MPD arc thrusters. 
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